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Introduction 

One of the most widely used programs for transonic unsteady aerodynamic 
analysis is the LTRAN2 code of Ballhaus and Goorjian.l That code is used to 
solve the low frequency approximation of the transonic small disturbance (TSD) 
equation. Steady state boundary conditions are used at the airfoil, in the 
wake, and on the computational boundaries. 

Use of the low frequency approximation and steady state airfoil and wake 
conditions limit the frequency of unsteady motion that can be analyzed with 
LTRAN2. Houwink and van der Vooren 2 extended the range of applicability of 
LTRAN2 by adding unsteady terms to the airfoil and wake boundary conditions; 
the resulting code was termed LTRAN2-NLR. Hessenius and Goorjian 3 added a 
time derivative term in the downstream far-field condition as well as unsteady 
airfoil and wake conditions. Their code, LTRAN2-HI, has been validated in the 
transonic range by a series of comparisons with experimental data. 

Although adding unsteady terms to the airfoil and wake boundary conditions 
extended the range of applicability of LTRAN2, use of the low frequency 
approximation of the TSD equation still limits its application to relatively 
low frequency motions. The programs described in references 1-3 use Murman- 
Cole (M-C) type dependent spatial differencing^ that admits nonphysical 
expansion shock waves as part of the computed solutions. Using steady far- 
field conditions causes disturbances incident on the boundaries to be reflected 
back into the computational domain. This necessitated placing the computa- 
tional boundaries far enough from the airfoil that reflected disturbances did 
not reach the airfoil during the calculations. Having to place the boundaries 
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far from the airfoil increases the cost of using the programs because the flow 
field has to be computed at an increased number of grid points. 

To remove some of the limitations described above, a new code, XTRAN2L, 
has been developed at the NASA Langley Research Center. It was developed by 
modifying LTRAN2-NLR. The M-C differencing was replaced with Engquist-Osher 
(E-0) monotone differencing. 5-5 E-0 differencing does not admit expansion 
shocks as part of the solution and increases code efficiency by allowing larger 
time steps to be used in the time-marching solution. The low frequency 
limitation was removed by adding the capability of solving the complete TSD 
equation. The final modification that is discussed in the present work is the 
implementation of nonreflecting far-field boundary conditions that are 
consistent with the complete equation. 

Edwards et al.? added the capability of including aeroelastic effects in 
the time-marching XTRAM2L solutions, and Seidel et al .8 made an extensive 
study to determine optimum methods for distributing cartesian grids. The 
details of those efforts may be found in the cited material. 

Unsteady Transonic Small Disturbance Equations 
The codes described in references 1-3 are used to solve the low frequency 
approximation of the transonic small disturbance (TSD) equation 

xt = ^xx + fyy (1) 

where <j> is a disturbance velocity potential normalized by cU62/3 } c i s 
airfoil chord, 6 is airfoil thickness ratio, and U is freestream speed. The 
spatial coordinates, x and y, and time, t, are normalized by c, c/6^/ 5 , and 
uj-l, respectively, where u> is the frequency of unsteady motion. The 
coefficient A = 2kMoo2/<s2/3 w ^ ere i s free-stream Mach number, and the 


reduced frequency k = wc/U. In references 1 and 3, 

B = (1 - M»2)/62/3 _ M oo ni( y + i )«j> x , where the choice of the exponent 
m is arbitrary. Ballhaus and Goorjian made m a function of such that the 
critical pressure coefficient, C p *, predicted by (1) matched the exact 
isentropic C p *. Hessenius and Goorjian used m = 2 (Spreiter scaling). In 
reference 2, B = (1 - NU 2 )/^ 2 /^ - M«> 2 (y* + lHx> where 
Y* = 2 - (2 - y)M«» 2 

Solution Algorithm 

In the codes described in references 1-3, solutions of (1) are obtained 
using the alternating-direction-implicit (ADI) scheme described in reference 
9. Solutions are advanced from the nth level in time to level n+1 using the 
following two-step procedure 
x-sweep: 


6 U. . - d) 1 ? .) = D f. . + 6 d) 1 ? . 

At x' v i,j m,o x l.o yy*i,j 

y-sv/eep: 

6 - (J). .) = is ( - 4>? . ) 

At X IS 1,J y i,j' 2 yy'M,j y i,j 


where $ is an intermediate level potential. From references 1 and 9, 
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(LTRAN2-NLR) 
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The mixed difference operator, D x , is constructed to maintain conservation 
form. Murman-Cole (M-C) spatial differencing used in LTRAN2 , LTRAN2-NLR, and 
LTRAN2-HI results in the following form for D x fjj: 
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(3b) 


It has been shown that M-C differencing allows stable, entropy-violating 
expansion shocks to be computed as part of the numerical solution. 

Reference 6 also showed that M-C differencing can trigger numerical 
instabilities that cause large errors in the calculated aerodynamic loads. 

Such a case was calculated using LTRAN2-NLR for flow over an MBB-A3 airfoil 
oscillating in pitch about its leading edge at = 0.8, k = 0.2. A time 
step of kAt = 1° was used, and the pitching motion was defined by an unsteady 
angle of attack a(t) = - 0.5° + 0.5°sin(kt). Figure la shows that the 
steady flow field is mixed subsonic/supersonic with a shock wave of moderate 
strength located at approximately 65 percent chord. Figures lb-lf show that 
during the airfoil oscillation an instability is triggered at the lower leading 
edge that causes the calculations to diverge. When the monotone differencing 
scheme of Engquist and Osher 5 is used, expansion shocks are not admitted as 
part of the computed solution, and the calculations remain stable when methods 
using M-C differencing have begun to diverge. 



Engquist-Osher Differencing 


The Engquist-Osher (E-O) scheme was first used in implicit algorithms by 
Goorjian and Van Buskirk^ who tested the method using a modified LTRAN2 
code. Similar modifications were made to LTRAN2-NLR at the NASA Langley 
Research Center. To incorporate the E-0 method into the ADI procedure required 
the following differencing during the x-sweep: 


At 5 x^i ,j " * 1 # j) “Vi -1/2, j + 6 yy*i,J 

where 
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The complete set of difference equations using the monotone differencing for 
the x-sweep are presented in Appendix A. 

To demonstrate the effect of the monotone differencing on numerical 
stability, the case of the oscillating MBB A-3 airfoil was recalculated. The 
pressure distributions. Figure 2, show that the numerical solution remained 
stable for the duration of the calculations. Goorjian and Van Buskirk reported 
that for some cases, they were able to increase kAt (and hence code efficiency) 
by factors of up to 10 and still maintain stability. 

Algorithm for the Complete TSD Equation 
The ability to treat unsteady motions of all frequencies was obtained by 
adding the capability to solve the complete TSD equation 

C* tt + Mxt = B<i>xx + byy 

where (5) 

r _ k 2 M 2 

C '? 7J 


Solutions for are advanced from time level n to time level n+1 using the 

following ADI method of Rizzetta and Chin** 

x-sweep: 


tv D f. .+ 6 . 

At x' i , j Y i,r x i , j yy y i,.i 


(6a) 


or 


~ 6 (<)). . - <j>’} .) = U f. . + S <b? . 
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(6b) 


y- sweep: 

_C 
At' 
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For the x-sweep, the algorithm is the same as that for the low frequency 
equation. Including <j>^ also requires an extra level of computer storage — 



levels n+1 , n, and n-1 versus levels n+1 and n. The difference approximation 
for the complete TSD equation are presented in Appendix B. Since the x-sweep 
is unchanged, only the difference equations for the y-sweep are presented. 

Nonreflecting Boundary Conditions 

The steady state far-field boundary conditions used in LTRAN2 cause 
disturbances incident on the grid boundaries to be reflected back into the 
computational domain. Thus, the boundaries had to be placed far away such that 
reflected disturbances did not reach the airfoil during the calculations and 
cause errors in the computed solution. Kwak 12 incorporated the nonreflecting 
far-field boundary conditions of Engquist and Majda^ into LTRAN2 which 
allowed a reduction in the physical extent of the computational grid and saved 
between 10 and 24 percent in computer time. The boundary conditions of 
reference 13 are not compatible with (5). Nonreflecting far-field boundary 
conditions that are consistent with the complete TSD equation are presented in 
this section. 

Assuming B to be locally constant, the transformations 




(7) 

(8) 
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In untransformed coordinates, (8) becomes 

F F + D)<f> t + F*x + f^y + If = 0 (9) 

Allowing x to approach - °° in (9) with y finite, the following first order 
plane wave condition at the upstream boundary was obtained 

14 + TF^t " *x = 0 (10) 

Similarly, letting x + 00 with y finite resulted in the downstream condition 


h~ F + TF^t + *x = 0 


( 11 ) 


As y ■*■ + “> with x finite, the following conditions at the top and bottom 
boundaries were obtained 

§<l> t ±<i> y =0 (12) 

where + and - represent the top and bottom boundaries, respectively. Using 
= f(r-t), — a solution of (7) that represents outgoing waves — to replace 


<(>t b y - qy (7 - F^~ 1<f> x’ ( 12 ) became 

7T - S = 0 (13) 

The boundary conditions in (13) were used in all numerical experiments. The 
difference equations for (10), (11), and (13) are presented in Appendix C. 

When C = 0, ( 10) - ( 13) reduce to the boundary conditions for the low frequency 
equation. 12 The nonreflecting boundary conditions are summarized in 
Figure 3. 

One test of the boundary conditions was in the calculation of unsteady 
forces on an NACA 64A010 airfoil pitching harmonically (about its quarter 
chord) ±0.25 degrees (°) about a 0° mean angle of attack at M„ = 0.825 
and k = 0.5. For the steady flow, an embedded shock wave is located at 



approximately 62 percent chord. A reference solution was calculated for four 
cycles of oscillation (360 steps per cycle) on a 113 x 97 grid (in x,y) that 
extended 200c in x and 709c in y. The grid was reduced to 88 x 65 
(-3.8c £ x £ 3.5c, y £ 9.3c), and the calculations were made first using 
steady-state far-field boundary conditions and then using (10), (11), and (13) 
at the boundaries. As shown in Figure 4, when the steady conditions were used, 
disturbances reflected from the boundaries caused the calculated lift to 
deviate significantly from the large grid value. When the nonreflecting 
boundary conditions were implemented, most of the waves incident on the 
boundaries were absorbed, and the small grid results showed good agreement with 
the reference calculation. Those results are also shown in Figure 4. Compared 
with the time required to generate the large grid solution (3215 seconds on a 
CDC CYBER 173), using the new boundary conditions on the small grid resulted in 
a 44 percent savings in computer time (the small grid solution required 1815 
seconds) . 

A second test was to calculate the unsteady force response for a flat 
plate airfoil with a pulse in angle of attack a. The calculations were made 
for M«> = 0.85 on an 80 x 61 grid that extended ±20c in x and ±25c in y. 

Using the pulse/transfer function technique described in Reference 8, the 
frequency response function for the unsteady lift curve slope c £q was 
calculated with and without the nonreflecting boundary conditions. In the 
pulse/transfer function technique, after a was increased smoothly and rapidly 
to a maximum and returned to its initial value, calculation of the unsteady 
forces were continued until those forces returned to their starting values. A 
Fast Fourier Transform (FFT) of the lift coefficient c £ was then divided by 
the a FFT to obtain the frequency response function for c £ q . A flat plate 
airfoil was used to allow comparisons of the forces calculated using XTRAN2L 
with those predicted using the exact kernel function method of Bland. 15 
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Figure 5a shows a comparison of the unsteady forces calculated using steady 
state conditions on the computational boundaries with the forces obtained using 
Bland's kernel function method. Below k * 0.5, the finite difference results 
have spurious oscillations due to disturbances reflected from the boundaries. 
When the nonreflecting boundary conditions were used (Figure 5b), the reflected 
disturbances were small, and good agreement with the exact solution was 
obtained. 

Concluding Remarks 

A new computer program, XTRAN2L, for transonic unsteady aerodynamic 
analysis has been developed at the NASA Langley Research Center. It is a 
modification of the LTRAN2-NLR code. The monotone differencing method of 
Engquist and Osher was used to replace the Murman-Cole type dependent 
differencing scheme. That resulted in a code that is more robust and 
more efficient, and the new differencing method does not admit nonphysical 
expansion shocks as part of computed solutions. The capability of analyzing 
airfoils undergoing motions of all frequencies was obtained by adding a general 
frequency term to the transonic small disturbance (TSD) equation. Solutions of 
the complete TSD equation are advanced through time using the alternating- 
direction-implicit method of Rizzetta and Chin. Nonreflecting boundary 
conditions that reduced disturbances reflected from the computational 
boundaries were implemented. This allowed the boundaries to be moved closer to 
the airfoil and thus further increased program efficiency. 
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APPENDIX A 


DIFFERENCE EQUATION FOR THE ENGQUIST-OSHER METHOD 


The difference equation that results when Engquist-Osher (E-O) monotone 
differencing is used in the x-sweep of the solution procedure is presented in 
this Appendix. When the E-0 method is used, the finite difference 
approximation of (4) becomes 
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where 
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The difference equation for the y-sweep is unchanged 



APPENDIX B 


DIFFERENCE EQUATIONS FOR THE COMPLETE TRANSONIC SMALL DISTURBANCE EQUATION 


The difference equations for the complete transonic small disturbance 
(TSD) equation are discussed in this Appendix. Since the numerical procedure 
for the x-sweep is the same for the low frequency and the complete TSD 
equations, only the difference equation for the y-sweep (6c) is presented 
here. That equation is 
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The tri diagonal coefficients are the same as for the low frequency equation 
with the exception of the term added to B. and the (2d. 1 ? . - 


term added to Dj. 
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APPENDIX C 


DIFFERENCE EQUATIONS FOR THE NONREFLECTING BOUNDARY CONDITIONS 

The difference equations for the nonreflecting boundary conditions are 
presented in this Appendix. The upstream and downstream boundary conditions 
are implemented during the x-sweep of the ADI procedure. They are applied 
midway between the extreme and adjacent grid columns at a time level halfway 

between level n and the intermediate level at which <j> is defined (n+1/2). At 
the upstream boundary, i = 1, the plane wave condition is 
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At the downstream boundary, i = IMAX (the maximum streamwise grid 
location), the nonreflecting condition is 
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E i c - Cl/2, j - Cl,J> 

At the lower boundary, j = 1, the boundary condition that is imposed is 
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At the top boundary, j = JMAX, 
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B i -1/2, j-1/2 


1 - M 


- - . + 4 ,? . , - 4>? n 

~ 2 VH i,j M-1,J M,J-1 v i-l,j-l) 


The tri diagonal coefficients are 


A _ 1-1/2, j-1/2 1 

j " x i - x i-i " y j - >j-i 

= a i -1/2, j-1/2 + 1_ 


J x i' x i-l y j“ y j-l 


C. = 0 
J 


v- 


* n+ ; . + . . <{, n+ ; . . 

y i-i,j y i-i,j y i,j y i-i,j-i 


- * 


i-l.j-l 


* n • i 

i ,J-1 


y j • y j-l 


a. 


i -1/2, j-.l/.2 ( ^n + + n ,"+} _ n+1 n Ji } 

x i - x^ v M,j-l v i,j v i -1 , J -1 -1 , J v i-l,j-l v i-l,J ; 


At the point (i,j) = (2,1) 


(4 - a. . (<j> )" + VJ = 0 

^y ; i,J+l i»J+l vq x ; i,j+l 


(Cll) 


where <{> x and <J>y were approximated with backward differences, and in the 
definition of ai,j+i 
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B, 


1 - M 


.n 


i,j+l * 2/3 ’ • • -/ x . . , - x. 


The difference equation for (Cll) is 

,n+l n+1 n n 

*j,in - *i..i * *i..hi ~ *1 , j 


■- H f (T * + 

i+1 " x i -1 


.n+1 


.n+1 


,n 


a** ' * .»» • * , , m ,n 

a. . . * *1 .3*1 ' Vi, j+l 


y j+l - y j 


i J+l 


x i " x i-l 


(C12) 


Using the upstream boundary condition 


^t ^i -1/2, j+l ' b i -1/2. j+l f *x -1/1, J+l = 0 


n+1/2 


the relationship 
,n+l 


‘i’i*} i +1 = 4>i i +1 + (1 + — -1/2, j+l — j-l/j b i-l/2,j+l At . , n n+1 . 

1 1,J+ 1 ,>J+1 x i- x 1-l ’ 1 x i - *,_! ,( *i -1, j+l ' 

was substituted in (C12). The tridiagonal coefficients at (i ,j) = ( 2 ,1) then 

became 

A. = 0 

J 


B, = 1 

J y j+l - y j 


•j -5 — u + (1 + *i.-l/2,J+l at r l (1 . t Vl/2,j+lV | 

J y j+l y j x i - x i-l x i - X i-1 1 U X. - X._ x >1 

1 " X i-1 x i ~ x i-l j x i - x i _! ^i -1 J+l 


a. . , 

D = 

j x 


,n 


f. . . - d>v . 

.J+l 9 1..1 

y j+i - y j 

At (>,j) - (2.JMAX), 




(C13) 
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where <|> x and <fy were approximated with backward and forward differences, 
respectively, and in a-jj_i 


B. . , 
i ,J-1 


i - m; 
T27T 


M f(y* + 


1) 


A 0 A 0 

♦i+1,3-1 *i-l,J-l 


i+1 


‘ x i-l 


Combining (C13) and the upstream condition 


(» t ) 


n+1/2 
i —1/2 ,j—l 


+ b 


i-l/2J-l 




n+1/2 

i-l/2,j-l 


The tri diagonal coefficients became 


A. = i ♦ -. Vfc l- [1 ♦ (1 ♦ W.J-l“ r l ( i 

3 - >4-i x i - x i-i 1 x i - x i-i 

i 

B, = 


b i-l/2,j-l At 
x i - x i-l 


>] 


J y j ' y M 

C. =0 


D - n . fl . b i-l/2,j-l At r l {1 b i-l/2J-l 

" Vi-i 1 x i - x i-i 


.• i A t 


n 

1-1J-1 


A 11 A n 

*i,J I *U-1 

y j - ^-1 
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